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O Three dimensional studies of convection in deep spherical shells have been used to test the hypothesis 

that the strong jet streams on Jupiter, Saturn, Uranus, and Neptune result from convection throughout the 

^ molecular envelopes. Due to computational limitations, these simulations must be performed at parameter 

1—1 settings far from Jovian values and generally adopt heat fluxes 5-10 orders of magnitude larger than the 
planetary values. Several numerical investigations have identified trends for how the mean jet speed varies 
with heat flux and viscosity, but no previous theories have been advanced to explain these trends. Here, 
we show using simple arguments that if convective release of potential energy pumps the jets and viscosity 
T-H damps them, the mean jet speeds split into two regimes. When the convection is weakly nonlinear, the 

equilibrated jet speeds should scale approximately with F/f, where F is the convective heat flux and v 
is the viscosity. When the convection is strongly nonlinear, the jet speeds are faster and should scale 
approximately as {F/vY^"^. We demonstrate how this regime shift can naturally result from a shift in the 
behavior of the jet-pumping efflciency with heat flux and viscosity. Moreover, both Boussinesq and anelastic 
simulations hint at the existence of a third regime where, at sufficiently high heat fluxes or sufficiently small 

. . viscosities, the jet speed becomes independent of the viscosity. We show based on mixing-length estimates 

^ ^ that if such a regime exists, mean jet speeds should scale as heat flux to the 1/4 power. Our scalings provide 
a good match to the mean jet speeds obtained in previous Boussinesq and anelastic, three-dimensional 

\^ simulations of convection within giant planets over a broad range of parameters. Whcin extrapolated to 

^ the real heat fluxes, these scalings suggest that the mass-weighted jet speeds in the molecular envelopes 
of the giant planets are much weaker — by an order of magnitude or more — than the speeds measured at 
cloud level. 
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At the cloud levels near ^Ibar pressure, numerous 
east- west (zonal) jet streams dominate the meteorology of 
the giant planets Jupiter, Saturn, Uranus, and Neptune, 
but the depth to which these jets extend into the interior 
remains unknown. Endpoint theoretical scenarios range 
from weather-layer models where the jets arc confinc^d to 
a layer several scale heights deep to models where the 
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jets extend throughout the molecular envelope ('--^lO^km 
thick) on cylinders parallel to the rotation axis (for a 
review see Vasavada and Showman 2005). For Jupiter, 
in-situ observations by the Galileo probe at 7°N latitude 
show that the equatorial jet extends to at least ~20 bars 
(~150 km below the visible clouds) (Atkinson et al. 1997), 
and indirect inferences suggest that the jets at other lati- 
tudes extend to at least ~5-10 bars pressure (e.g., Bowl- 
ing 1995; Legarreta and Sanchez-Lavega 2008; Morales- 
Juberias and Bowling 2005; Sanchez-Lavega et al. 2008). 
At Neptune, gravity data suggest that the fast jets are 
confined to the outermost few percent of the planet's mass 
(Hubbard et al. 1991). Comparable data are currently 
lacking for Jupiter and Saturn but will be obtained by 
NASA's Juno and Cassini missions, respectively, in com- 
ing years. 

Three-dimensional (3B) numerical simulations of con- 
vection in rotating spherical shells have been performed 
by several groups to investigate the possibility that the 
jets on the giant planets result from convection in the 
interior. Both free-slip and no-slip momentum boundary 
conditions at the inner and outer boundaries have been 
explored; of these, the free-slip case — which allows the 
development of strong zonal flows — is most relevant to gi- 
ant planets. So far, such studies have neglected the high 
electrical conductivity and coupling to magnetic fields ex- 
pected to occur in the deep (^ 1 Mbar) planetary interior. 

In this line of inquiry, most studies to date make the 

Boussinesq assumption in which the basic-state density, 
thermal expansivity, and other background properties are 
assumed constant with planetary radius; the convection 
is driven by a constant temperature difference imposed 
between the bottom hot boundary and top cold bound- 
ary (Aurnou et al. 2008; Aurnou and Heimpel 2004; Au- 
rnou and Olson 2001; Christensen 2001, 2002; Heimpel 
and Aurnou 2007; Heimpel et al. 2005). These studies 
show that convection with free-slip boundaries can drive 
multiple jets with speeds greatly exceeding the convective 
speeds. Studies using thick shells tend to produce ~ 3-5 
jets (Am-nou and Olson 2001; Christensen 2001, 2002). 
When the shell thickness is only ^10% of the planetary 
radius, then at least under some parameter combinations, 
^15— 20 jets can occur, similar to the number observed on 
Jupiter and Saturn (Heimpel and Aurnou 2007; Heimpel 
et al. 2005). Nevertheless, many factors in addition to 
shell thickness can affect the number of jets. 

In real giant planets, the density and thermal expan- 
sivity each vary by several orders of magnitude from 
the cloud layer to the deep interior, and a new gener- 
ation of convection models is emerging to account for 
this strong radial variation in basic-state properties. Us- 
ing the anelastic approximation, which accounts for this 
layering, Evonuk and Glatzmaier (2006, 2007), Evonuk 



(2008) and Glatzmaier et al. (2009) present idealized 
two-dimensional (2B) simulations in the equatorial plane 
exploring hypothetical basic-state density profiles, with 
density varying by up to a factor of 55 across the con- 
vection zone. In contrast, Jones and Kuzanyan (2009) 
present 3B simulations using an idealized basic-state den- 
sity structure, with density varying by a factor of up to 
148, while Kaspi et al. (2009, 2010) present 3B simula- 
tions with a realistic Jovian interior structure, with den- 
sity varying by nearly a factor of 10** from the deep in- 
terior to the 1-bar level. These anelastic studies likewise 
suggest that the jets could penetrate deeply through the 
molecular envelope. 

A challenge with all of the above-described simula- 
tions is that, for computational reasons, they must be 
performed using heat fluxes and viscosities that differ 
greatly from those on Jupiter, Saturn, Uranus, and Nep- 
tune (Fig. 1). Thus, while simulations can produce jets 
with speeds similar to the observed values of ^100- 
200 m sec" ^ for some combinations of parameters (e.g., 
Aurnou et al. 2007, 2008; Aurnou and Olson 2001; Chris- 
tensen 2001, 2002; Heimpel and Aurnou 2007; Heim- 
pel et al. 2005; Jones and Kuzanyan 2009; Kaspi et al. 
2009), this does not imply that convection in the in- 
terior of real giant planets would necessarily produce 
jets with such speeds. In fact, depending on the pa- 
rameter combinations, simulations with free-slip bound- 
ary conditions can produce jets that equilibrate to mean 
speeds^ ranging over many orders of magnitude, from ar- 
bitrarily small (less than Imsec"^) to lOOOmsec"^ or 
more. Assessing the likely wind speeds in the molecu- 
lar envelopes of the real giant planets — and determining 
whether their observed jets can be pumped by convection 
in their interiors — requires the development of a theory 
that can be extrapolated from the simulation regime to 
the planetary regime. 

Currently, however, there is no published theory that 
can explain the jet speeds obtained in simulations with 
free-slip boundaries nor their dependence on heat flux, 
viscosity, and other parameters. Several investigations 
have presented scaling laws describing how the mean 
zonal-wind speeds vary with control parameters when 
no-slip boundary conditions are used, as potentially rel- 
evant to Earth's outer core (e.g., Aubert 2005; Aurnou 
et al. 2003). These arc essentially theories for the mag- 
nitude of wind shear in the fluid interior for cases where 
the zonal velocity is pinned to zero at the boundaries. 
However, these scaling laws are not applicable to the gi- 
ant planets, where the fluid at the outer boundary can 
move freely, lacks a frictional Ekman layer, and exhibits 



^Calculated in a suitably defined way, such as yj2K/M, where 
K and M arc the total kinetic energy and mass; this is essentially 
the mass- weighted characteristic wind speed of the fluid. 
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strong jets. Several attempts have also been made to 
quantify how the convective velocities scale with param- 
eters, but — regardless of the boundary conditions — these 
relationships cannot be applied to the jet speeds because 
the convective and jet speeds can differ greatly and their 
ratio may depend on the heat flux and other parameters. 

The goal of this paper, therefore, is to develop scal- 
ing laws for how the jet speeds depend on heat flux and 
viscosity that explain the simulated behavior within the 
simulated regime and, ideally, allow an extrapolation to 
real planets. The simulations themselves make a number 
of simplifications (e.g., ignoring magnetohydrodynamics 
in the deep interior at pressures exceeding ~1 Mbar), but 
in our view building a theory of this idealized case is a 
prerequisite for understanding more realistic systems. 

We first quantify the degree of overforcing in current 
studies, since this issue has received little attention in the 
literature (Section 2). Next, we quantify the dependence 
of the convective speeds on planetary parameters and 
compare them to results from an anelastic general circu- 
lation model from Kaspi et al. (2009) (Section 3). Armed 
with this information, we construct simple scalings for 
the characteristic jet speeds in three regimes. In the first 
two regimes (Section 4), the viscosity is large enough so 
that viscous damping of the jets provides the dominant 
kinetic-energy loss mechanism. Christensen (2002) sug- 
gested the existence of a third regime where the jet speeds 
become independent of the viscosity; we construct possi- 
ble scalings for this regime in Section 5. In Section 6, we 
combine the three regimes and discuss extrapolations to 
Jupiter. Section 7 concludes. 

2. Degree of overforcing 

For numerical reasons, current 3D simulations of con- 
vection in the giant planets must use viscosities many or- 
ders of magnitude larger than the molecular viscosities. 
This results from the coarse grid resolution in the models: 
to be numerically converged, such a model must have con- 
vective boundary layer and convective plume thicknesses 
of at least a gridpoint, and this requires very large viscosi- 
ties. Given the enhanced damping of the kinetic energy 
implied by this constraint, such simulations can achieve 
wind speeds similar to those on the giant planets only by 
adopting heat fluxes orders of magnitude too large. 

Moreover, even in the absence of such a viscous damp- 
ing, enhanced heat fliixes are computationally necessary 
simply to achieve spin up within available integration 
times. To illustrate, imagine a giant planet that initially 
has no winds, and ask how long it would take to spin 
the zonal jets up to full strength if the jets penetrate 
through the molecular region (down to ^1 Mbar pressure 
on Jupiter and Saturn). This corresponds to a kinetic 
energy per area of ~ u^Ap/g, where u is the characteris- 
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Fig. 1. Illustration of the parameter space associated with 
convection in rotating spherical shells. The modified-flux 
Rayleigh number Ra*p (abscissa) and Ekman number E (or- 
dinate) can be viewed as non-dimensionalized heat flux and 
viscosity, respectively. The Prandtl number, giving the ratio 
of viscosity to thermal diffusivity, constitutes a third dimen- 
sion. Published simulations access a small region of parameter 
space with values of Rap and E typically 10~® or greater, but 
the actual values for Jupiter are ~10~^'' or less. It is unknown 
whether convection at Jupiter-like values of Rap and E would 
produce Jupiter-like wind speeds. 



tic jet speed (~ 30 m see" ^ on Jupiter), g is gravity, and 
Ap ~ 1 Mbar is the pressure thickness across which the 
jets are assumed to extend. On giant planets, the work 
to spin up the winds comes from the convection, driven 
by internal heat fluxes ranging from 0.3 Wm~^ or less for 
Uranus and Neptune to 5 W m^^ for Jupiter. Given these 
small available fluxes, the characteristic pumping time of 
the jets is ~10^ years. This is a lower limit, because it as- 
sumes that the efficiency of converting the heat flux into 
kinetic energy is close to 100%; in reality, the efficiency 
will be less, and much of the power produced will be used 
to resist frictional losses in the simulations. In contrast, 
published state-of-the-art high-resolution numerical sim- 
ulations, if expressed in dimensional units, typically inte- 
grate months to years (e.g., Aurnou et al. 2008; Aurnou 
and Olson 2001; Ghristensen 2001, 2002; Heimpel and 
Aurnou 2007; Heimpel et al. 2005; Jones and Kuzanyan 
2009; Kaspi et al. 2009). Thus, the heat fluxes must be 
enhanced by at least 3-5 orders of magnitude simply to 
allow the simulations to spin up within achievable inte- 
gration times. 

Most convection papers cast the equations in nondi- 
mensional form, rendering unclear the actual degree to 
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which these sinmlations are overforced. One can quantify 
the degree of overforcing as follows. In these simulations, 
the specified control parameters are the Rayleigh number 
Ra, the Ekman number E, and the Prandtl number Pr, 
defined as follows: 



Ra = 



agATD^ 



E = 



Pr=- (1) 



where a, g, k, and u are the thermal expansivity, grav- 
ity, thermal diffusivity, and kinematic viscosity. D is the 
thickness of the convecting layer, AT is the imposed tem- 
perature drop across this layer, and fl is the planetary 

rotation rate;. 

The heat &ux, which is an output, is typically cast 
as a Nusselt number Nu, which gives the ratio of the 
total (convective plus conductive) heat flux to a reference 
conductive flux that would be conducted across a system 
with the same temperature drop and thickness: 



Nu = 



J tot 
^cond 



DF,, 



KpCp 



AT' 



(2) 



where Ftot is the total flux, the thermal conductivity is 
KpCp, and Fcond = KpCpAT/ D. At Rayleigh numbers far 
exceeding the critical value, the total heat flux approxi- 
mately equals the convective heat flux. 

So how may we determine the dimensional heat flux 
for a nondimensional simulation with specified Ra, E, 
Pr, and Nu! Using (1) and (2) shows that 



-Ftot — 



pcp E^RaNu 



ag 



Pj,2 



(3) 



Because convection is driven by buoyancy associated 
with a heat flux, an alternate Rayleigh number, deflned in 
terms of the heat flux rather than a temperature contrast, 
is useful: 

agF^otD^ 



Rap 



pCpUK'' 



(4) 



Christensen (2002) argued that, in the limit of small 
diffusivities, the fluid behavior should become indepen- 
dent of the diffusivities, and he defined a modified heat- 
flux Rayleigh number that is independent of these diffu- 
sivities. This can be constructed by multiplying Eq. (4) 
by E'^Pr^'^, giving^ 



Ra% = 



agFt, 



pCp 



(5) 



^Christensen (2002) made alternate definitions of Ra^ and Nu 
that included a non-dimensional factor of rj (the ratio of the inner- 

to-outcr radius adopted in the simulations) in the numerator of 
Eqs. 2 and 4. We forgo this factor here. 



In the limit of small diffusivities, one thus expects that 
the behavior should scale with Ra*p rather than Ra. Cast- 
ing Eq. (5) as an expression for heat flux, we have 



pCp 



ag 



Ra*p. 



(6) 



For a giant planet, all the dimensional parameters in 
Eq. (6) are known. Consider Jovian values of f7 = 1.74 x 



10- 



sec 



D 



2 X 10"* km, Cp = 1.3 x lO^Jkg-iK-i, 



and g = 26msec~^. For interpreting Boussinesq simula- 
tions in the context of giant planets, an ambiguity is that 
the thermal expansivity and density vary greatly from 
atmospheric values at the top to liquid values in the inte- 
rior. Using atmospheric values corresponding to a pres- 
sure of ~100bars, where the temperature is ~600K, we 
have p ~ 4kgm~^ and a ~ 1.6 x 10~^K~^. Consider- 
ing the simulation in Christensen (2001), where Pr = 1, 
E = 3x 10"^, Ra = 10^, and Nu « 10, we then have heat 
fluxes of ~10^ W m~^. However, given that the vast bulk 
of the interior behaves as a liquid, much more appropriate 
values are p ~ 10^ kgm~^ and a ~ 10~^ K~^, relevant to 
the bulk of Jupiter's interior. These values imply fluxes 
of ^10*^ IQ-'^^Wm^^. Similar estimates imply that the 
simulations of Heimpel and Aurnou (2007) and Aurnou 
et al. (2008) adopt a heat flux of ~10^ Wm'^. Aurnou 
et al. (2007) 's simulations of Uranus and Neptune adopt 
a heat flux of ~10^ Wm'^. 

Given the great variation of density and thermal ex- 
pansivity with radius on a giant planet, future studies 
will increasingly adopt models that incorporate the full 
radial variation of these parameters. This will require 
extending the definition of Rayleigh number. We adopt 
the following definition, in which a mass- weighted average 
is performed over the radially varying quantities (Kaspi 
et al. 2009): 

Ral = J^('^), (7) 



where the mass- weighted average is deflned as (...) = 
J{...)pr'^ dr/ J pr'^ dr. 

Anelastic simulations that adopt realistic radial pro- 
files for the density and thermal expansivity lack the am- 
biguity described above. Kaspi et al. (2009) force the flow 
not through a constant-temperature boundary condition 
but by introducing internal heating and cooling. Integrat- 
ing their heating profile over the depth of the planet gives 
an effective heat flux with a maximum of ~10^Wm~^ 
in the deep interior with values decaying to zero in the 
atmosphere. Similarly, Jones and Kuzanyan (2009) es- 
timate a dimensional heat flux of ~10®Wm~^ for their 
simulations, broadly consistent with the estimates made 
here. 

The implication is that the simulations in Christensen 
(2001, 2002), Aurnou and Olson (2001), Heimpel and Au- 
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rnou (2007), Aurnou et al. (2008), Kaspi et al. (2009) and 
related studies are overforced by ~6-10 orders of mag- 
nitude. The overforcing compensates for the numerical 
need to use high diffusivities and achieve a steady state 
over reasonable time scales. In non-dimensional terms, 
the numerical models adopt Ra*p of lO"'^ or greater while 
Jupiter's Ra*p value is ~10" 



-14 



3. Convective velocities and buoyancies 

In a giant planet, convection in the interior transports 
the interior heat flux. To order-of-magnitude, one thus 
expects that 

F - pCpw5T, (8) 

where F is the convected heat flux, w is the characteris- 
tic magnitude of the vertical convective velocity, and 6T 
is the characteristic magnitude of the temperature dif- 
ference between a convective plume and the surrounding 
fluid. 

In a rapidly rotating, low-viscosity fluid, the scaling 
for convective velocities is typically written (e.g., Boub- 
nov and Golitsyn 1990; Fernando et al. 1991; Golitsyn 
1980, 1981; Stevenson 1979)^ 

It is known that Eq. (9) provides a reasonable represen- 
tation of convective velocities for laboratory experiments 
in rotating tanks (e.g., Boubnov and Golitsyn 1990; Fer- 
nando et al. 1991; Golitsyn 1981). Several authors have 
also suggested that Eq. (9) is relevant for the dynamo- 
generating regions of planetary interiors, where a three- 
way balance between buoyancy, Coriolis, and Lorentz 
forces is often assumed (e.g., Starchenko and Jones 2002; 
Stevenson 2003). To our knowledge, however, no broad 
assessment of its accuracy has been made for convection 
in the molecular interior of a giant planet. To do so, 
we performed three-dimensional numerical simulations of 
convection in Jupiter's interior using an anelastic gen- 
eral circulation model based on the MITgcm (Kaspi et al. 
2009). The radial dependence of gravity, basic-state den- 
sity, compressibility, and thermal expansivity correspond 
to a realistic Jupiter interior structure calculated with 
the SCVH equation of state (Saumon et al. 1995). Sim- 
ulations with a broad range of parameters were explored 
(Table 1). All of our simulations are in a rapidly rotating, 
low- viscosity regime, with geostrophic balance holding on 
large scales (see, e.g., Kaspi et al. 2009, Fig. 6). 

^Equation (9) can be heuristically derived by assuming that the 
convective buoyancy force per mass, ga ST, is balanced by the verti- 
cal Coriolis force due to the turbulent eddy motions, Qu' , where u' 
is the horizontal eddy velocity. Assuming that the turbulent eddy 
motions are approximately isotropic («' ^ w) and invoking Eq. (8) 
then leads to Eq. (9). 
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Fig. 2. Top: Mass-weighted convective vertical wind speed 
(calculated as the deviation of vertical wind speed from its 
zonal mean) versus modified heat-flux Rayleigh number for a 
range of anelastic simulations. The velocities are expressed 
as a convective Rossby number, defined as velocity divided 
by 57_D. Dashed line is the scaling of Eq. (10) with 7 = 1. 
See Table 1 for definitions of the symbols. Bottom: Radial 
variation of the horizontally averaged velocity components for 
an anelastic simulation with parameters Ra'p — 2.89 x 10~^, 
i5 = 1.5 X 10~*, and Pr = 10. Depicts deviation of vertical 
velocity from its zonal mean (dashed), deviation of the zonal 
velocity from its zonal mean (dashed-dotted), the square root 
of their zonally averaged correlation [u'w'Y^'^ (dotted), and 
the scaling from Eq. (9) evaluated using the radially varying 
thermal expansivity a, basic-state density p, heat flux F, and 
gravity. Note that all the anelastic simulations in this paper 
use the full Jovian interior structure with a factor of ~10* 
variation in density from top to bottom. For more details on 
the experiments see Kaspi (2008) and Kaspi et al. (2009). 
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10 
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1.07, 1.82,3.15,6.8,6.33, 
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square 


blaclc 
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square 
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triangle 
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left-pointing 
triangle 


black 
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circle 
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0.8,1,2,3,4, 
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diamond 


red 
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8,8.5,10, 12 


10 
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34.9,27.9,23.9,18.6, 17, 
16.2,14.3,13^,11.4, 127, 
11.2,10.1,9.26,9.39,7.58, 
10.6 
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10 
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Table 1. List of numerical integrations, using the model of Kaspi et al. (2009), that are presented in this paper. Each row lists 

a sequence of simulations. Columns 2, 3, and 4 give the control parameters, and column 5 gives the global-mean, mass- weighted 
Rossby number. The last two columns give the symbol type and color used to present each sequence of simulations in Figs. 2, 
5, and 8. All the anelastic model integrations in this paper adopt the full Jovian radial interior structure of gravity, density, 
thermal expansivity, compressibility, and other thermodynamic properties from Kaspi et al. (2009), corresponding to a factor 
of ~10'* variation of density from top to bottom. 
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Our simulations show that Eq. (9) provides a good 
approximation over a wide range of heat fluxes and rota- 
tion rates. Figure 2 shows the mass- weighted mean con- 
vective velocities from the Kaspi ct al. (2009) anclastic 
simulations (presented as a convective Rossby number, 
namely convective vertical velocity divided by flD) ver- 
sus modificd-flux Raylcigh number and compares it to 
Eq. (9). For the vast majority of simulations, the sim- 
ulations differ less than a factor of ~3 from the scaling, 
indicating that Eq. (9) is quite accurate. This analysis 
complements Kaspi et al. (2009) 's finding that, within a 
specific simulation, Eq. (9) captures the approximate ra- 
dial dependence of the convective velocity across the in- 
terior, over which a/p varies by a factor of ~ 10^ (Fig. 2, 
bottom). Together, these results suggest that Eq. (9) is 
valid not only for the mass-weighted mean velocities but 
locally within the fluid. 

Laboratory experiments suggest that a prefactor should 
exist on the right side of Eq. (9) with a value close to 2 
(Boubnov and Golitsyn 1990; Fernando et al. 1991; Golit- 
syn 1981); however, our simulations are better matched 
when no prefactor is used. 

We can nondimensionalize the convective velocities us- 
ing a convective Rossby number, Roconv = w/ {Q,D). The 
convective flux, F, equals F^ot minus the conductive flux. 
The conductive flux is only important near the critical 
Rayleigh number, and in this limit it becomes i^cond from 
Section 2. Making the approximation that the conductive 
flux equals i^cond, Eq. (9) nondimensionalizes to yield 

Roconv ~ liRa*F - Ra*p^''y/^, (10) 

where Ra*p^^* is the critical modified flux Rayleigh num- 
ber for the onset of convection. Thus, at the criti- 
cal Rayleigh number, the expression properly predicts 
zero convective velocities, and at greatly supercritical 
Rayleigh numbers it predicts Rocom ~ 

7(i?a^)i/2. We 

have introduced a dimensionless prefactor 7 for general- 
ity, but consistent with Fig. 2, we will assume 7=1 when 
performing numerical estimates. 

Like Eq. (10), many of our subsequent scalings will 
depend on the difference between the actual and critical 
Rayleigh numbers. For notational brevity, we therefore 
define this difference as 



Rcltp = RaX 



(11) 



Equations (8-9) imply that the fractional density as- 
sociated with convective plumes should scale as 

a^T4^V'^ (12) 

Aubert et al. (2001) proposed an alternate scaling for 
the convective velocities in the limit of negligible vis- 
cosities, Roconv ~ {Ra*p)^^^. They performed rotating 



laboratory experiments in liquid gallium and water, and 
showed that this expression provides a reasonable fit to 
the convective velocities inferred for their experiments. 
In the numerical simulations of Christensen (2002) and 
Christensen and Aubert (2006), the poloidal component 
of the velocity field (which includes the convective veloc- 
ities) depends on the Ekman number, but at small Ek- 
man number seems to be converging toward an asymp- 
totic dependence that is reasonably well represented by 
this 2/5 scaling (with a prefactor of 0.5). In our case, 
Eqs. (9)-(10) provides a slightly better fit to the convec- 
tive velocities, although the 2/5 scaling (with a prefactor) 
is also adequate. Conversely, overplotting Eqs. (9) (10) 
against Christensen (2002) 's simulation data shows that 
they match essentially as well as the 2/5 scaling. 



4. Scaling for the jet speeds 

a. Experimental data 

Our goal is to understand the physical processes gov- 
erning the global-mean jet speeds as a function of heat 
fiux and viscosity (i.e., Rossby number as a function of 
Ekman and Rayleigh numbers) for low-viscosity, rapidly 
rotating convection in spherical shells. To characterize 
this functional dependence requires numerous (~100 or 
more) numerical integrations so that the available param- 
eter space is adequately sampled. While it is numerically 
possible in three-dimensional simulations to reach Ekman 
numbers as low as 3-4 x 10~^ (e.g., Aurnou et al. 2008; 
Heimpel and Aurnou 2007; Heimpel et al. 2005; Jones and 
Kuzanyan 2009), this requires intensive computation and 
precludes the exhaustive survey of parameter space re- 
quired to develop scaling laws. Instead, our strategy is to 
adopt more modest Ekman numbers (extending to 10~^ 
to 10~^), which ensure that the dynamics are still in the 
low-viscosity, rapidly rotating regime (with geostrophic 
balance holding on large scales) yet allow numerous sim- 
ulations to be performed. Here we present over 100 sim- 
ulations. In addition, we extensively use the dataset of 
Christensen (2002), who systematically characterized the 
global-mean Rossby numbers as a function of Ekman and 
Raylcigh numbers in Boussinesq simulations with Ekman 
numbers as low as 10~^. Despite being far from Jovian 
parameter values, all the simulations we present are in the 
appropriate rapidly rotating, low-viscosity regime, with 
geostrophic balance holding at large scales. 

These numerical experiments demonstrate that the 
mean jet speeds depend significantly on the control pa- 
rameters Ra*p and E. Figure 3 depicts the global-mean, 
mass- weighted Rossby numberi?o = U/ {^D), where D is 
the thickness of the model domain and U is the domain- 
mean, mass-weighted wind speed, as a function of the 
control parameters in the Boussinesq simulations from 
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Fig. 3. Dependence of mass-weighted, domain-averaged jet 
speed on parameters from published simulations. Depicts the 
Rossby number (color and contours) as a function of modi- 
fied flux Rayleigh number (abscissa) and Ekman number (or- 
dinate) for Boussinesq simulations from Christensen (2002) 
(top) and anelastic simulations as described in Table 1 using 
the model of Kaspi et al. (2009) (bottom). Note that con- 
tours are evenly spaced in log(_Ro) and that all figures in this 
paper use identical contour values and colorbars, facilitating 
intercomparison. In both panels, the plusses denote the loca- 
tions of individual simulations. The range of Ekman numbers 
plotted is different in the two cases; anelastic cases are com- 
putationally more demanding and hence were run at larger 
Ekman numbers than the Boussinesq cases. The small-scale 
structure results from interpolation onto a fine grid of the 
coarsely spaced experimental points and may not be robust. 
The robust feature is the overall trend of positive slopes of 
constant- iio contours, with mean slopes close to one, and Ro 
values ranging from ~0.I on the right to ~0.000f on the left. 
Note that, in both sets of simulations, the contours are widely 
spaced at Rossby numbers exceeding ~0.02 (right side of plot) 
but become tightly spaced at Rossby numbers less than ~0.02 
(left side of plot), suggesting a regime transition. 



Christensen (2002) (top panel) and the anelastic simu- 
lations described in Table 1 using the model of Kaspi 
et al. (2009) (bottom panel). The domain-averaged jet 
speeds range across three orders of magnitude even within 
the limited parameter space explored in these models, 
with some parameter combinations producing Jupiter- 
like speeds and others not. Rossby numbers range from 
^0.0001 to 0.1, corresponding to domain-mean wind 
speeds of ~1 to lOOOmsec"^. In both cases, larger Ra*p 
and smaller E promote faster wind speeds. This makes 
sense qualitatively because, for a given rotation rate and 
planetary size, larger Ra*p implies larger heat flux (i.e., 
stronger forcing of the flow), while smaller E implies 
smaller viscosity (i.e., weaker damping of the flow). Kaspi 
et al. (2009) found that constant-wind-speed contours ex- 
hibit a slope (in the \og{Ra*p)-log{E) diagram) of ap- 
proximately 5/4. The data from Christensen (2002) plot- 
ted in Fig. 3 likewise indicate that the wind-speed con- 
tours in his case exhibit a mean slope of approximately 
1. For both sets of simulations, Fig. 3 also shows that the 
contours tend to be widely spaced toward the right and 
tightly spaced toward the left, with an approximate tran- 
sition at Ro « 0.02. This suggests a transition between 
two regimes. 

The processes that determine how mean jet speeds de- 
pend on Ra*p and E are not understood, however, and 
so the trends in Fig. 3 remain unexplained. Develop- 
ing such an understanding is important because rotat- 
ing spherical-shell convection is an inherently interesting 
physics problem, and also because extrapolation into the 
jovian parameter regime can only be performed once a 
theory for the Ra*p and i?-dependence of the mean jet 
speeds has been developed. In the absence of such a the- 
ory, it is unknown whether convection at Jupiter-like val- 
ues of Ra*p and E would generate Jupiter-like jet speeds. 
Explaining Fig. 3 is the core goal of this paper. 

b. Regime I: Strongly nonlinear regime 

Here, we construct a simple scaling theory based on en- 
ergetic arguments to attempt an explanation for the jet 
speeds in the regime of fast jets {Ro ^ 0.02 in Fig. 3, cor- 
responding to jet speeds exceeding ^ lOOmsec^^), which 
we call Regime I. 

First consider the forcing. A convecting fluid parcel 
traveling at a vertical speed w releases potential energy 
per mass at a rate P ~ gwSp/p, where p is density and Sp 
is the magnitude of the density contrast between plumes 
and the background environment. Using the equation 
of state Sp = apST, where p represents a background 
density and 5T is the magnitude of the temperature con- 
trast between plumes and the background environment, 
we can write P ~ gwa ST. This potential energy is pri- 
marily converted into convective kinetic energy. Some 
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fraction e of this energy is used to pump the zonal jets; 
we provisionally assume e is constant but return to this 
assumption in Section 4 c. Integrating over the planetary 
mass, the total power (in W) pumped into the jets is then 
approximately 



-ftot ~ 47re / aSTgwpr^ dr, 



(13) 



where the mass of an infinitesimal spherical shell of radial 
thickness dr is Airpr^ dr. Both ST and w are a priori 
unknown, but they are related by the fact that convection 
transports the planet's interior heat flux (Eq. 8), which 
implies that w 5T ~ F/(pCp). Thus, the power available 
to pump the jets can be written 



Pu 



47re 



agF 



■ dr. 



(14) 



Energy loss occurs through friction, which we assume 
acts as a diffusive damping of the winds with a viscosity 
V (for the simulations, this would be the model viscos- 
ity that enters the definition of the Rayleigh and Ekman 
number; see Eq. 1). The power per mass dissipated by 
viscosity is approximately v'S/'^v? ^ vk^u^, where k is the 
wavenumber of the structures dominating the dissipation, 
that is, the dominant wavenumber of V^u^. This implies 
a total rate of kinetic energy loss given by 



7 2 2 2 I 

fk u pr dr. 



(15) 



Now we equate energy gain (Eq. 14) and loss (Eq. 15) 
to obtain an expression for mean jet speed. In the case of 
a Boussinesq fluid, we note that a and p are constants and 
that F , g, Cp, and u are approximately constant (within 
a factor of ^-^3) and can be pulled out of the integral, 
yielding"* 

\ PCpV J 



(16) 



The case of an anelastic fluid is complicated by the large 
variation of a and p with radius, and even u may vary 
significantly from the top to the bottom of the domain 
(Kaspi et al. 2009). However, in the Kaspi et al. (2009) 
simulations, the mean jet speeds are relatively constant 
(within a factor of ^ 3) throughout most of the domain; 
only the top few % of the fiuid mass experiences signifi- 
cant wind shear along the Taylor columns. To obtain a 
crude expression for the mass-weighted mean jet speed 
in the molecular region, we therefore assume that u is 
constant, allowing us to write: 



-'(;) 



1/2 



dr 



J pr"^ dr 



1/2 



(17) 



similar equation was derived, albeit with different arguments, 
by IngersoU and Pollard (1982, Eq. 20). 
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Fig. 4. Prediction of the scaling law for Regime I given 
in Eq. (18), in which the jet-pumping efficiency e is assumed 
constant. Depicts contours of Rossby number versus modified 
flux Rayleigh number (abscissa) and Ekman number (ordi- 
nate). Adopts kD — Stt and e = 0.3 (see text). Note the 
similarity of the slopes and values of the contours with the 
simulation results plotted in Fig. 3, particularly for Rossby 
numbers greater than ~0.02. 

Let us nondimensionalize these expressions so that 
they may be compared with the simulation results shown 
in Fig. 3. For both cases, we use Ro = u/{VtD) and the 
definition of E given in Eq. (1). Using the appropriate 
definition of the modified flux Rayleigh number yields the 
same nondimensional expression for both the Boussinesq 
and anelastic cases. Again approximating the convective 



flux F as the total flux minus Fr, 



Ro: 



eV2 
kD 



Ra^ 

E 



, we obtain 

1/2 



(18) 



where Ra'^ 



Ra% 



Ra*^''^ (see Eq. 11). In Eq. (18), 



Ra*p is given by Eq. (5) for the Boussinesq case and by 
Eq. (7) for the anelastic case. Note that, in the highly 
supercritical regime in which we expect Eq. (18) to be 
valid, the critical Rayleigh number is generally negligible, 
so that Ra^ approximately equals Ra*p. 

Thus, this simple theory predicts that the characteris- 
tic jet speed is proportional to (F/z^)*/^, or equivalently 
that the characteristic Rossby number is approximately 
proportional to {Ra*p/Ey^^ . Increases or decreases in 
Ra*p and E by the same factor leave Ro unchanged; in 
other words, this theory predicts that constant-i?o con- 
tours should have a slope of one in the Ra*p-E plane. 
This explains the fact that the empirically determined 
slopes in the simulations are close to one (see Fig. 3). In 
our theory, this behavior occurs because the forcing de- 
pends on the flux in the same way as damping depends 
on viscosity (namely, linearly). Thus, equal increases (or 
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decreases) in the flux and the viscosity cause alterations 
to the forcing and damping that cancel out, leading to a 
mean jet speed that is unchanged. 

To further compare the theory to the simulations, 
Fig. 4 plots our solution (Eq. 18) as contours of Rossby 
number versus Ra*p and E. A choice of kD is necessary 
to evaluate Eq. (18). Kaspi et al. (2009)'s simulations 
show that, although the jets themselves are broad (with 
widths comparable to the domain thickness), the curva- 
ture of the wind V^u, and thus the dissipation, exhibits a 
dominant wavelength of about D/5 (see their Fig. lOd). 
High-wavenumber structure in the flow shear also occurs 
in Christensen (2002) 's simulations. To account for this 
high-wavenumber structure, we adopt kD « Stt, and we 
calculate Fig. 4 using this value and e = 0.3. 

Figure 4 reiterates the approximate agreements in the 
slopes of constant-i?o contours (compare Fig. 4 with 
Fig. 3). Moreover, the Ro values themselves also agree 
well, particularly in the regime of fast zonal jets where 
Ro ;> 0.02. Throughout this part of the parameter space, 
the predicted Rossby number (Fig. 4) matches those ob- 
tained in the simulations (Fig. 3) to within a factor of 
^2 at a given Ra*p and E, althoiigh the discrepancy ap- 
proaches an order of magnitude toward the left side of the 
plot (at Ro ^ 0.02) where the simulated Rossby num- 
bers seem to undergo a regime shift that is not captured 
by Eq. (18). Such a shift could result for example from 
a dependence of e or fc on Rap and E, a possibility not 
considered up to now. We return to this point in the next 
two subsections. 

c. Regime II: Weakly nonlinear regime 

For sufficiently small values of Ra*p, both the Boussi- 
nesq and anelastic simulations shown in Fig. 3 exhibit 
a regime shift where the dependence of Rossby number 

on Ra*p and E becomes steeper than the {Ra*p/Ey^'^ 
dependence discussed in the previous subsection. As 
Fig. 3 shows, at Rossby numbers smaller than ~0.01- 
0.02, the constant-i?o contours becomes closely spaced, 
with a spacing indicating that a scaling {Ra*p/E)^ with 
^ sa 1 might provide an approximate fit. There is also 
an indication that the slopes of the constant-i?o contours 
decrease toward the left side of the plot, particularly for 
Christensen's simulations. 

These stronger dependences of Rossby number on Ra*p 
suggest a regime shift where different processes set the 
characteristic Rossby number at low Rossby number than 
at high Rossby number. In the strongly nonlinear regime 
explored in the previous subsection, the Rayleigh num- 
ber is strongly supercritical, and the convection is chaotic 
and nonlinear. At lower Rayleigh numbers, when the 
Rayleigh number is only modestly greater than the crit- 
ical value for convection, the convection is laminar and 



more spatially organized. Here, we explore how this tran- 
sition in convective behavior might lead to the regime 
shift in the jet speeds seen in Fig. 3. 

To do so, we consider an alternate approach based the 
zonal momentum balance. Consider a cylindrical coor- 
dinate system whose axis aligns with the rotation axis, 
with coordinates (s. A, z) corresponding to cylindrical ra- 
dius (i.e., distance from the rotation axis), longitude, and 
distance above or below the equatorial plane, respectively. 
The zonal-mean zonal momentum equation then reads (cf 
Kaspi et al. 2009) 

^ + 2f2t^-|- • (ptZv) -I- • {pu'w') = vV'^u (19) 

where u is the zonal wind, Vg is the cylindrically radial 
wind component (that is, the velocity component away 
from the rotation axis), overbars denote zonal means, 
primes deviations therefrom, and p is the mean density 
(constant in the Boussinesq case and a specified function 
of radius in the anelastic case). On the left side, the sec- 
ond term is the Coriolis acceleration. The third term 
represents advection due to the mean-meridional fiow, 
and the fourth term represents the acceleration caused 
by eddies (i.e., Reynolds stress convergences). Generally 
speaking, for the geostrophic jovian regime, the Reynolds 
stress term dominates over the advection term (e.g. Kaspi 
et al. 2009). 

We now average the equation in z (i.e. along the direc- 
tion of Taylor columns) . The Coriolis accelerations cancel 
out, because mass continuity prohibits any net, column- 
averaged motion toward or away from the rotation axis. 
Coupled with the free-slip boundary conditions and sym- 
metry of the zonal-wind structure about the equatorial 
plane, the ^-averaging also removes the ^-components 
of the divergence and Laplacian terms. Neglecting the 
mean-flow advection terms, the column- averaged, steady- 
state balance becomes 



s^p ds 



(s u'v'J) 



d_ 
'ds 



Id. 

S OS 



(20) 



which states that, in steady state, viscous drag associated 
with shear of the mean zonal wind du/ds balances jet 
acceleration caused by eddies. Balances analogous to this 
expression have been written, for example, by Busse and 
Hood (1982), Busse (1983a), Busse (1983b), and Cardin 
and Olson (1994).^ 

We here \isc this balance to achieve a simple expres- 
sion for the jet speeds. Approximate the left-hand side 
as ku'v'g, where k is the wavenumber (in cylindrical ra- 
dius) over which u'v'g varies significantly. Furthermore, 



^An analogous balance was considered by Aubert et al. (2001) 
and Aubert (2005) but for the no-slip boundary condition where 
boundary-layer friction plays the dominant role in the damping. 
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suppose that both u' and v'^ scale like our expected con- 
vective velocity scale, so that u'v'^ ~ Cw"^, where w is the 
convective velocity (e.g., from Eq. 9) and C is a correla- 
tion coefficient equal to one when u' and v'^ are perfectly 
correlated (i.e. when eastward u' always occurs with out- 
ward v'g and vice versa) and equal to zero when u' and v'^ 
exhibit no correlation. We approximate the right-hand 
side as vk^u. The momentum balance then becomes 



C 



w 
vk 



(21) 



Nondimensionalizing, the Rossby number associated with 
the jets is 

Ro^^"^ (22) 

An equation analogous to this was derived by Christensen 
(2002, his Eq. 3.12). Inserting our expression for the con- 
vective velocities (Eq. 10) yields 



Ro' 



C7^ Ra^ 
~kD^' 



(23) 



Thus, this theory predicts that, if the degree of cor- 
relation between u' and v'^ is independent of the control 
parameters (i.e., if C is constant), then the domain- mean 
jet speed scales as the convective flux over the viscos- 
ity, or equivalently the mean Rossby number scales as 
{Ra*p — iiajf "* )/E. Away from the critical Rayleigh num- 
ber, increases or decreases in Rap and E by the same fac- 
tor leave Ro unchanged; therefore, as with the theory pre- 
sented in the previous subsection, Eq. (23) predicts that 
constant-i?o contours should have slopes close to one in 
the logarithmic Ra*p E plane. Near the critical Rayleigh 
number, however, the predicted slopes differ from one. 
If the critical modified flux Rayleigh number depends on 

to a positive power less than one, then the constant- 
Ro contours increase in slope, becoming more vertical; if, 
however, Ra*"^^ depends on £ to a power greater than 
one, then the constant- i?o contours decrease in slope, be- 
coming more horizontal. 

To evaluate Eq. (23) quantitatively, we require an 
expression for iJaJf'*. By performing many simula- 
tions, Christensen (2002) determined empirically the crit- 
ical Rayleigh number for each Ekman number he ex- 
plored.^ Cardin and Olson (1994) performed a linear 
instability analysis of rotating convection in a spheri- 
cal shell and found that, at ^ 1, the critical value 



® Christensen found that the critical values of a modified 
Rayleigh number Ra* = RaE'^Pr-^ are 0.001005, 0.002413, 
0.006510, and 0.016790 for Ekman numbers of 10"^, 3 x 10-^, 
10~*, and 3 x 10~*, respectively. This is equivalent to modified 
flux critical Rayleigh numbers, Ra*^"*, of 1.005 x 10"®, 7.2 x 10"*, 
6.5 x 10-''', and 5.04 x 10"^ for those same Ekman numbers, respec- 
tively. Equation (24) gives values that agree with these to within 
~10%. 



of the ordinary Rayleigh numbcir is a constant times 
£.-i7/i5p^4/3(;L Pr)-''/3. Using the relationships be- 
tween Ra and Ra*p (see Section 2), and noting that the 
Nusselt number equals one at the onset of convection, 
shows that the critical modified fiux Rayleigh number 
should then scale as E'^^/'^^Pr-'^/^il + Pr)-*/^. Most 
of Christensen's simulations adopt a Prandtl number of 
one and we find that his empirically determined critical 
Rayleigh numbers are well matched by the expression'' 



Ra*p"'^ « 20i;28/i5. 



(24) 



Given this dependence, contours of constant Rossby num- 
ber in the logarithmic Ra*p-E plane should decrease in 
slope as one approaches the critical Rayleigh number (i.e., 
as one moves toward lower Rap values). This is quali- 
tatively consistent with the behavior exhibited by Chris- 
tensen's simulations (see Fig. 3). 

d. Combining Regimes I and II 

We have discussed two distinct regimes: (i) a regime 
of fast zonal winds and strongly supercritical convection 
where the simulated jet speeds seem well explained by 
the assumption that e is constant (Regime I), and (ii) a 
regime of slow zonal winds and weakly supercritical con- 
vection where the jet speeds are approximately explained 
by the assumption that the correlation between the u' 
and v'g velocity components is strong and roughly con- 
stant (Regime II). These led to two distinct scalings for 
the jet speeds: Eqs. (18) and (23) for the two regimes, 
respectively. Two issues now arise: First, how do we 
combine these regimes? Second, why should energetics 
and momentum considerations lead to different scalings? 
Both approaches should, in principle, yield the same an- 
swer within any given regime. 

The resolution to both issues lies in the dependence 
of the correlation coefficient C and the convective jet- 
pumping efficiency e on the control parameters. Con- 
sider, for example, the correlation coefficient. When the 
Rayleigh number is weakly supercritical and the flow is 
geostrophic, the convection forms broad, laminar convec- 
tion rolls, which in a sphere become tilted in the prograde 
direction (Busse 1983a, b, 2002; Busse and Hood 1982; 
Christensen 2002; Kaspi 2008; Sun et al. 1993; Tilgner 
and Busse 1997, and others). This leads to highly cor- 
related velocity components u' and v'^, as illustrated in 
Fig. 5 {top left): ascending fluid parcels move prograde 



'^Dormy et al. (2004) found theoretically that the critical 
Rayleigh number should scale as E~*/^ at low Ekman number, 
which would imply that iJo^''* should scale as E^/^. However, 
this scaling does not fit Christensen's empirically determined criti- 
cal Rayleigh numbers as well as Eq. (24), so we retain scaling (24) 
for the purposes of this paper. The choice has a negligible influence 
on our results. 
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Fig. 5. Top: Snapshots in the equatorial plane of two 
anelastic, 3D simulations from Kaspi (2008) illustrating how 
the correlation between the convective velocity components 
depends on the supercriticality. Each slice shows half of 
the equatorial plane of one simulation (note however that 
the simulations each span 360° of longitude). Color depicts 
streamfunction and arrows denote velocity component in the 
equatorial plane. The left simulation is weakly supercritical; 
the strong correlation between outward and eastward veloc- 
ity components is obvious. The right simulation is strongly 
supercritical; the correlation between the outward and east- 
ward velocity components is weaker because of the complex, 
turbulent convective structure. Bottom: Depicts correlation 
coefflcient C versus Rap for a range of anelastic simulations*. 
As expected qualitatively from the top panels, C decreases 
with Rap. Symbols are defined in Table 1. 



while descending parcels move retrograde. As a result, 
at weakly supercritical Rayleigh numbers, the correla- 
tion coefficient C should be close to one, at least outside 
the tangent cylinder. On the other hand, when the con- 
vection is strongly supercritical, the convective structure 
is chaotically time-dependent and spatially complex; the 
velocity components u' and v'^ are positively correlated 
in some regions but negatively correlated in others. This 
is illustrated in Fig. 5 {top right). The net correlation is 
usually still positive, but the large degree of cancellation 
implies that, at strongly supercritical Rayleigh numbers, 
the correlation coefhcient should drop significantly below 
one (Christensen 2002). As a result, one expects the cor- 
relation coefficient C to vary slowly with Ra*p near the 
critical Rayleigh number but to decrease more rapidly 
with increasing Ra*p at sufficiently supercritical Rayleigh 
numbers. This behavior can be seen in correlation coef- 
ficients calculated for a range of simulations versus Ra*p 
(Fig. 5, bottom). 

Next consider the jet-pumping efficiency e. The scaling 
(18), wherein the Rossby number scales with {Ra*p/E)^/^ , 
assumes that e is constant. While this assumption seems 
to work well at sufficiently supercritical Rayleigh num- 
bers, it must fail in the weakly supercritical regime 
where m' and are highly correlated. This is because, 
in this weakly nonlinear regime, the mean jet speeds 
scale quadratically with the convective velocities, and the 
power exerted to pump the jets scales as the Reynolds 
stresses times the jet speeds, namely as the fourth power 
of the convective velocities. The dependence of this 
quantity on the heat flux differs from the dependence of 
the convective potential-energy release on the heat flux. 
Since e is the ratio of these quantities, e cannot be con- 
stant in this regime. 

A regime shift analogous to that seen for the jet speeds 
in Fig. 3 — where Rossby numbers scale approximately 
with Ra*p / E at low Rossby number and approximately as 
{Ra*p/E)^/^ high Rossby number — can occur if the cor- 
relation coefficient depends weakly on Ra*p at low Ra*p 
yet strongly on Ra*p at high Ra*p. Likewise, it can occur 
if the jet-pumping efficiency depends strongly on Ra*p at 
low Ra*p yet weakly at high Ra*p. 

To quantify these arguments, let us relate e and C to 
each other. By definition, the efficiency e is the ratio of 
work done to pump the jets to the work made available 
by convection: 



m4V • (pu'v') 
gwa ST 



(25) 



where it is understood that both the numerator and the 
denominator represent global averages. Expressing the 
numerator as Ckw^, the denominator as gaF/pCp (using 
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Eq. 8), and non-dimensionalizing, we obtain 
RoCkDRol^^^ 

€ ~ A . 



(26) 



Using Eq. (22) for the Rossby number yields the expres- 
sion 



.Rot 



ERa^ ' 



(27) 



If we adopt Roc 
obtain finally 



j{Rapy^^ for concreteness, we 



^2 iR<4 

' E 



(28) 



Therefore, e and C cannot simultaneously be constant 
when Ra*p or E are varied. Holding one constant requires 
the other to become a function of Ra*p and E. 

Let us now make the simplest possible assumption re- 
garding this regime shift — we postulate a regime with 
C ~ constant = 1 at low Rayleigh number and a regime 
with e « constant = Cmax at high Rayleigh number. 
Quantitatively, there is no rigorous expectation that C 
or e need be constant, nor (as described previously) is 
this assumption actually necessary for a regime shift to 
occur. Nevertheless, C and e have upper limits of 1, so if 
some process causes them to increase with increasing (or 
decreasing) Ra*p, they might naturally plateau — at least 
over some range of parameter space — upon approaching 
their upper limits. Still, future theoretical work on what 
sets the Ra*p- and i?-dependence of the correlation coef- 
ficient and jet-pumping efficiency is warranted. 

Given these postulates, Eq. (28) implies that 



C = 
and 
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max7 E; 

-f-^E 



(29) 



(30) 



To smoothly transition between the regimes, we adopt 
the expression 



C 



7 

tmax 



Raj 



A\ 1/2 



n -1 



(31) 



which is merely the inverse of one over the value of C from 
the first regime plus one over the value of C from the sec- 
ond regime. This gives behavior equivalent to Eq. (29) in 
the two limits with a smooth transition in between. Given 
this expression, e can then be evaluated from Eq. (28), 
yielding 




10"^ 10 ^ 10 ' 10 



1.000 



0-100 



o 



o. 



0.010 - 



0.001 




lO"-" 10 



Fig. 6. Top: The correlation coefficient C, the ratio be- 
tween the Reynolds stress u'v's and the product of the root- 
mean-square velocity magnitudes. Bottom: e, the fraction of 
potential energy released by convection that is used to pump 
the jets. Blue triangles depict values we calculated rigorously 
from anelastic simulations at i5 = 4 x 10"*, while red squares 
denote values Christensen (2002) calculated rigorously from 
Boussinesq simulations at E = 10~^. In both cases, values 
represent flow behavior outside the tangent cylinder. Curves 
plot Eqs. (31)-(32) at an Ekman number of 4 x lO"* (solid 
blue) and 10~^ (dashed red) assuming 7 = 1 and tmax = 0.3, 
with Ra*p'^'^ evaluated using 

2Q£;28/i5 fgj. ^jjg Boussinesq sim- 
ulations and 0.4i5^*/'^''' for the anelastic simulations (the dif- 
ference resulting from the differing Prandtl numbers; see Ta- 
ble 1). 
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Figure 6 plots e and C versus Ra*p from available data 
and compares them to Eqs. (31)-(32). In particular, the 
blue triangles denote values of C and e that we rigorously 
computed^ from the three-dimensional convective veloc- 
ity and entropy fields for a sequence of anelastic simula- 
tions with E = A X 10~^. Additionally, the red circles 
denote values of C rigorously calculated by Christensen 
(2002) for a sequence of his simulations at E = 10~^ (see 
his Table 3; note that no data on e are available for Chris- 
tensen's simulations.) Blue and red solid curves denote 
the predictions of Eqs. (31)-(32) for i; = 4 x 10"^ and 
10~^, respectively, and can be directly compared to the 
symbols of the same color. 

Key points are as follows. First, the data indicate that, 
for both sets of simulations, C decreases with Ra*p, while 
e increases with Ra*p. Second, and more importantly, 
the slopes of the trends change across the range of Ra*p 
explored. The trends of C in both sets of simulations de- 
pend relatively weakly on Ra*p toward the left but seem 
to steepen toward the right. For example, Christensen's 
data indicate that a factor-of-two change in Ra*p cause a 
1.28-fold change in C at low Ra*p but a 2.7-fold change 
in C at high Ra*p. Conversely, the data suggest that e 
depends on Ra*p strongly at small Ra*p but more weakly 
at larger Ra*p. As Eq. (28) shows, these trends must 
go hand-in-hand; consistency between the two quanti- 
ties requires the trend in C to steepen if that in e flat- 
tens. Moreover, for both C and e, these changes in the 
slopes with Ra*p are precisely the ingredients that can 
generate a regime shift wherein Rossby numbers depend 
strongly on Ra*p at low Ra*p (as postulated for Regime II 
in Section 46) but more weakly on Ra*p at high Ra*p (as 
postulated for Regime I in Section 4a). Third, the data 
suggest that, at constant Ra*p, the correlation coefHcient 
increases with increasing E. 

The theoretical curves for C and e match the data 
surprisingly well (Fig. 6). All the qualitative trends dis- 
cussed above for the data are predicted by the theory: C 
decreases and e increases with Ra*p; the slope of C{Ra*p) 
steepens toward the right while that of e{Ra*p) flattens 
toward the right; and at constant Ra*p the values of C in- 
crease with E. The theory also predicts that at constant 
Ra*p, the values of e should decrease with increasing E, 



° The correlation coefficient and jet-pumping efficiency were cal- 
culated from the anelastic simulations as 



C ■- 



and 



[(«'2 t;'2)l/2] 



[«iV • (p«'v')] 



[gwa 5T] 

where [...] = J {...)ps ds/ J psds indicates a mass- weighted average 
over the equatorial plane. 




Fig. 7. Prediction of our scaling law, Eq. (33), that com- 
bines Regimes I and II. Depicts contours of Rossby number 
versus modified flux Rayleigh number (abscisssa) and Ekman 
number (ordinate). Adopts kD — 5n, 7=1, and tmax ~ 0.3 
(see text). Note the similarity of the slopes and values of the 
contours with the simulation results plotted in Fig. 3. Con- 
tours are evenly spaced in log(-Ro), with contour values iden- 
tical to those in Figs. 3 and 4. For reference, the dashed line 
depicts the critical modifled flux Rayleigh number given by 
Eq. (24). Above the dashed line, convection does not occur. 



which can be tested with future simulations. Moreover, 
for the parameters chosen (7 = 1 and e^ax — 0.3), the 
numerical values of the theoretical e and C curves agree 
relatively well with the simulation data, particularly for 
the anelastic simulations at £' = 4 x lO"**. On the other 
hand, the theoretical correlation coefficients at E ^ 10^^ 
overpredict those for the corresponding Boussinesq sim- 
ulations by up to a factor of several, particularly at the 
higher-_Ra|^, values. Of course, the technique adopted to 
smooth between the two regimes (cf Eqs. 31-32) will af- 
fect the degree to which the curves agree with the data, 
so the details of the comparisons should be considered 
tentative. 

What are the implications of this discussion for the 
jet speeds? Inserting our expression for e into Eq. (18) 
or our expression for C into Eq. (23) leads to a single, 
self-consistent expression for the Rossby number span- 
ning both regimes. For simplicity, we adopt Eqs. (29)- 
(30) rather than their smoothed counterparts, leading to 
the expression 



Ro — —— min 



2 RcL*p 



Ra*j^"^ 



.1/2 
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RcLp 



E 



1/2-1 



(33) 

which is equivalent to the minimum of Eqs. (18) and (23). 
Thus, arguments based on energetics and momentum ar- 
guments are now self-consistent: both approaches pre- 
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diet that the Rossby number seales as {Ra*p — Ra*p"*')/E 
at weakly supercritical Rayleigh numbers but as {Ra*p — 
Ra*p"^ /Ey^"^ at strongly supercritical Rayleigh numbers. 

Figure 7 plots the jet speed predicted by Eq. (33) as 
contours of Rossby number versus Ra*p and E, with a 
contour spacing that is even in log(i?o) and contour val- 
iica identical to those in Figs. 3 and 4. Regime I, in which 
the jet-pumping efficiency e is constant, lies toward the 
lower right, while Regime II, in which the correlation coef- 
ficient C is constant, lies toward the upper left. The tran- 
sition between the regimes manifests visually as a jump 
in the spacing of the constant-i?o contours and occurs 
at a Rossby number of ~0.02 for the parameters chosen, 
roughly consistent with the regime transition in the sim- 
ulated data (see Fig. 3). Within Regime II, the slopes of 
the contours decrease as the Rossby number is decreased 
and the modified flux Rayleigh number approaches the 
critical value. A similar decrease in the slope of the con- 
tours at small Rossby number is evident in Christcnsen 
(2002) 's data plotted in Fig. 3. Moreover, throughout 
most of the parameter space accessed by the simulations, 
the predicted Rossby number (Fig. 7) matches those ob- 
tained in the simulations (Fig. 3) within a factor of ~2. 

5. Scaling for the jet speeds in Regime III: the 
asymptotic regime 

The simulated behavior in Fig. 3 is reasonably well 
explained by the theory presented in Section 4, where 
Rossby numbers scale approximately with {Ra*p—Ra*p^^*') / E 
at low Rayleigh numbers (Regime II) and as {Ra*p/ E)^^"^ 
at high Rayleigh numbers (Regime I). However, Chris- 
tcnsen (2002) argued that at the highest Ra*p he ex- 
plored for each value of E (or cquivalcntly the lowest 
E values explored for each value of Ra*p) his simulations 
approached an asymptotic regime where the mean equi- 
librated jet speed became independent of the values of 
the diffusivities and depend only weakly on the heat flux. 
This behavior manifests in Fig. 3{top) as both a widen- 
ing of the constant-_Ro contours and significant increase 
in their slope beyond one toward the far right edge of 
the plot, neither of which is predicted by the scalings in 
Section 4. Based on an empirical fit to his simulation 
results, Christensen (2002) proposed that in this asymp- 
totic limit — which we call Regime III — the Rossby num- 
ber associated with the characteristic jet speeds scales 

9 



cally explained, we present here a scaling analysis based 
on mixing-length theory. We compare it to available sim- 
ulations and evaluate its implications for the Jovian pa- 
rameter regime. Nevertheless, we emphasize at the outset 
that the existence of the asymptotic regime is tentative, 
and further numerical work is required to confirm (or re- 
fute) its existence and determine its properties. This sec- 
tion is intended simply to offer some ideas on how such 
an asymptotic regime, if it exists, could work. 

We hypothesize that, if the viscosity is sufficiently 
small, the damping will be determined not by the molec- 
ular viscosity but by an eddy viscosity associated with 
turbulent diffusion. The concept of an eddy viscosity is 
relevant only when the convection is sufficiently supercrit- 
ical (so that the eddies behave in a nonlinear, turbulent 
manner), leading us to adopt the approach of Section 4a. 
To obtain an expression for mean wind speed, we again 
equate forcing (Eq. 14) and damping (Eq. 15) but use 
the eddy viscosity, feddy, in place of in Eq. (15). We 
adopt a standard mixing-length formulation for the eddy 
viscosity, i^eddy ~ 'f^^mbc' wherc w is the convective ve- 
locity (given by Eqs. 9-10) and fcmix is the wavenumber 
associated with the mixing length (e.g., a typical distance 
traversed by coherent convective plumes). In the case of 
a Boussinesq fluid, where a and p are constant and F, g, 
and Cp are approximately constant, this yields 



as 



Ro = 0.53{Rapy/^. 



(34) 



At present, however, a physical basis for this empirical 
fit is lacking. To see whether Eq. (34) can be theoreti- 



® Because he included an inner-to-outer radius ratio in his def- 
inition of Rap, Christcnsen (2002) quotes the relationship with a 
pre-factor of 0.65 rather than 0.53. 
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(35) 



where e is assumed approximately constant as in Section 

4a. For an anelastic fluid, where a and p vary greatly 
with radius, we follow the approach used in Section 4 a, 
yielding 



u ~ 



dr 



/(" 



^^^^\^^py\^ dr 



1/2 



(36) 



Let us nondimensionalize these expressions so that 
they can be compared with existing Boussinesq and 
anelastic numerical results and the empirical asymptotic 
scaling (Eq. 34) proposed by Christensen (2002). Defin- 
ing the Rossby number as Ro = u/{^D) and assuming 
that fcmix is constant and that Ra*p > Ra*"''^, Eqs. (35)- 
(36) can be nondimensionalized to yield 



Ro 



1/2 



{Rapy/^ 



(37) 



where Ra*p is given by Eq. (5) for the Boussinesq case 
and Eq. (7) for the anelastic case. F is a dimensionless 
constant that depends on the planet's basic-state radial 
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density structure and is given by 

r 1/2-1 
(/ pr^ drf/^ (j^r^ dr) 

r= , J/2 — (38) 

where /(r) is an order-unity dimensionless function that 
describes the radial dependence of the heat flux, F = 
Fo/(r), and Fq is a constant (independent of radius) 
that gives the characteristic value of the heat flux (e.g., 
halfway through the layer). Note that, for the case of 
constant p, a, g, flux, and Cp, then F = 1. Interestingly, 
for the Kaspi et al. (2009) model, T = 1.02 despite the 
large variation of a and p with radius.^" 

This simple theory therefore predicts that, in the 
asymptotic regime, the characteristic jet speed scales as 
F^/^, or equivalently the characteristic Rossby number 
scales as (i?a|^)i/*, with no explicit dependence on the 
molecular or numerical viscosity. Our exponent of 1/4 
compares favorably to — but is modestly steeper than — 
Christensen's empirically determined exponent of 1/5. 
Notably, our theoretical exponent is significantly smaller 
than that for Regimes I and II (Section 4), where, at 
constant E, jet speeds scale with F^/"^ or F depending 
on assumptions. In our asymptotic theory, the weak de- 
pendence on Ra*p occurs because the eddy viscosity de- 
creases with decreasing heat flux (unlike the behavior in 
Regimes I and II, where the viscosity is a fixed param- 
eter). As a result, a weakly forced fiow has extremely 
weak damping, whereas a strongly forced fiow has ex- 
tremely strong damping; this leads to an equilibrated jet 
speed that depends only weakly on heat fiux. Specifically, 
the kinetic-energy forcing scales with F, and because the 
convective velocities scale with F^/^ (Eq. 9), the eddy vis- 
cosity and therefore the kinetic-energy damping also scale 
with F^/^. This leads to an equilibrated kinetic energy 
scaling with F^/^ and a jet speed scaling with F^/"*. 

Figure 8 explicitly compares the theoretical scaling 
(Eq. 37) with the empirical scaling (Eq. 34) and the sim- 
ulation results of Christensen (2002) and Kaspi et al. 
(2009). The theoretical scaling is shown for a value of 
the prefactor e^/2(fcmixF')^/^7"^/^(fcF')"\ which we caU 
X, equal to 1. For this value, the magnitudes of the jet 
speeds predicted by our Eq. (37) match Christensen's fit 
within a factor of ~2 (compare solid and dashed lines in 
Fig. 8a). For the values of kD w Stt adopted previously, 
this would require A;„iix ~ (STr)^/!?, implying a rather 
short mixing length 27rA:~?^ of ^0.031?. This seems at 




10"" 10"' 10"" 10"^ 10"^ 

Ra*p 



Fig. 8. Comparison of our asymptotic scaling for the jet 
speeds, Eq. (37) (solid curves) with Christensen's empirical 
asymptotic fit from Eq. (34) (dashed curve) and simulation 
results from Christensen (2002) (panel a) and Kaspi et al. 
(2009) (panel b). Solid curves are evaluated using e = 0.3, 
7 = 1, and x ~ ^- See Table 1 for definitions of the symbols 
used in (6). 



'^'^Given the integrals in the numerator and denominator of F, 
one can show that V tends to be close to 1 for a wide range of 
possible radial dependences of the function agf /cp — even when this 
function varies in radius by orders of magnitude — as long as it has 
a smooth dependence, as occurs in the numerical models and in 
Jupiter. 
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least qualitatively consistent with the fact that the vor- 
ticity structure of convective plumes exhibits little coher- 
ence in cylindrical radius s (see Christensen 2002). The 
eddy diffusivitics implied by this choice of mixing length 
are similar to those suggested by previous authors, lend- 
ing some encouragement to this choice of prefactor.^^ 

Overall, the close agreement between the heat-flux de- 
pendences of our asymptotic scaling and Christensen's 
empirical fit is encom-aging and argues that the basic idea 
encapsulated by our scaling —that the damping depends 
on to a power modestly less than that of the forcing — is 
correct. 

At high Rayleigh numbers, both the Boussinesq and 

anelastic simulations (Fig. 8a and b, respectively) con- 
verge toward a trend similar to the shallow slope of the 
predicted asymptotic scaling. However, the simulations 
diverge from the asymptotic scalings at low Rayleigh 
numbers as the effects of the model diffusivities become 
strong. Christensen's fit, and the scaling presented in this 
section, are intended to describe the asymptotic behavior 
in the limit of negligible diffusivities (for discussion see 
Christensen 2002). 

When extrapolated to Jovian values of Ra*p, this 
asymptotic regime predicts that the mass-weighted mean 
velocities in the molecular envelopes of Jupiter and Sat- 
urn arc small, ^0.1-1 mscc"'^. We consider this issue 
further, and combine the three regimes into a single scal- 
ing, in the next subsection. 

6. Combining the three scalings 

We have derived scalings for the jet speeds in three 
regimes: two regimes in which the numerical viscosity 
dominates the damping (Regimes I and II in Sections 4 a 
and 4&, respectively) and another in which the viscosity is 
determined by turbulence, i.e., an eddy viscosity (Regime 
III in Section 5). Here, we combine the scalings. 

Our basic approach in deriving the scalings was to bal- 
ance forcing against damping; the damping represents 
a numerical (or molecular) viscosity in Regime I but 
an eddy viscosity in Regime HI. Here, we suppose that 
both molecular and eddy viscosities operate simultane- 
ously, and that whichever viscosity is larger will dom- 
inate. Because larger viscosity implies smaller equili- 
brated jet speeds, we therefore expect that the smaller 

^^The (dimensional) eddy diffusivity implied by our expression 
for vertical velocity (Eq. 9) and the condition x = 1 is I'eddy ^ 
eD(agF)^/^/[(pCj,f2)^/^7(fcZ))^]. For Jovian interior conditions 
and heat flux (a « IQ-^K"!, g « 26msec-2, F = 5Wm-2, 
pKi lOOOkgm-3, Cp « 1.3x10* J kg-iR-i, n = 1.74x 10-*sec-i 
D = 2 X 10* km), and e = 0.3 and kD = 5n as adopted previously, 
this implies values of ~ 20m^sec~^. This is roughly consistent 
with the value of ~10m^ sec~^ suggested by Starchenko and Jones 
(2002). 



of the two Rossby numbers will dominate. Combining 
Regimes I, II, and III, we can thus roughly say 
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(39) 



where we have assumed that k is constant and Cmax has 
the same values in the second and third terms on the 
right side. Neglecting Ra*p"^ in the expression for Regime 
I, the transition between Regimes I and HI occurs for 
Ekman numbers 



E, 



1/2 



tr 



^^krnixD 



(40) 



At a given Rap, Regime I (or II) occurs for large Ek- 
man numbers while the asymptotic regime occurs for 
small Ekman numbers, with a transition near E^J■. Anal- 
ogously, at constant E, Regime I (or II) dominates at 
small Ra*p while the asymptotic regime occurs at large 
Ra*p, with a transition near modified flux Rayleigh num- 
bers of rV^(fcmix-D)2i;2. 

To illustrate the combined scaling. Fig. 9 depicts 
Eq. (39) as contours of Rossby number versus Ra*p and 
E for the specific case Cmax = 0.3, 7=1, and kD = Stt 

(the same values as in all previous figures). The value 
of fciuix is chosen so that the prefactor in the asymptotic 

regime emax(fcmix£')^/^(7^''^fc^)~^ = X = 1- Regimes II 
and I lie in the upper left and exhibit sloped contours. 
The asymptotic regime lies at the lower right and ex- 
hibits vertical contours, consistent with the expectation 
that Rossby number is independent of E there. As ex- 
pected from Eq. (40), the transition between Regimes I 
and HI slopes gradually upward to the right; the asymp- 
totic regime extends to larger E when Ra*p is larger. For 
the value of fcmix adopted in Fig. 9, Christensen's simu- 
lations lie predominantly within Regimes I and II; how- 
ever, the predicted transition to the asymptotic regime 
occurs close to Christensen's highest-iJaJ. cases for each 
value of E. If Fig. 9 is correct, it would suggest that the 
asymptotic regime would become clear at Ekman num- 
bers an order of magnitude smaller than Christensen ex- 
plored for each value of Ra*p. Note, however, that in 
reality the transition between the regimes will probably 
occur smoothly across a broad strip straddling Etr{Rap); 
if so, it might be necessary to reach Ekman numbers as 
low as 10~^ to sec the asymptotic regime clearly. 

Now we examine Fig. 10a, which plots our combined 
scaling (Eq. 39) versus Ra*p for several values of E and 
compares them to Christensen (2002) 's Boussinesq sim- 
ulations. Curves and simulation data are shown for 
E = 10"^ (black), 3 x 10^^ (magenta), lO"'^ (blue), 
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Fig. 9. Predicted jet speeds for a scaling law that com- 
bines Regimes I, II, and III, given in Eq. (39), for the case 
Emax ~ 0.3, kD = Stt, and 7=1. The value of fcmix 
is chosen so that the prefactor in the asymptotic regime, 
X = ^lli^ikrai^DY^^ /{j^^'^kD), equals one. Depicts contours 
of Rossby number versus modified flux Rayleigh number (ab- 
scissa) and Ekman number (ordinate). In the Regimes I and 
II, contours tilt upward to the right, whereas in the asymp- 
totic regime they are vertical. For reference, the dashed line 
depicts the critical modified flux Rayleigh number given by 
Eq. (24). Above the dashed line, convection does not occur. 



Fig. 10. Predicted jet speeds for a scaling law that com- 
bines Regimes I, II, and III, given in Eq. (39), for the case 
Emax = 0.3, kD = Stt, and x = 1- Depicts Rossby number ver- 
sus modified fiux Rayleigh number for several values of E for 
our combined scaling (curves) and from simulations (symbols) . 
Panel a: Christensen (2002) 's Boussinesq simulations, at Ek- 
man numbers of 1 x 10"^ 3 x 10"^, 10"*, and 3 x 10"'', with 
scaling predictions evaluated at those same values (plotted in 
the corresponding colors). Panel h: Anelastic simulations at 
Ekman numbers of 1.5 x 10"*, 4 x 10"*, and 8 x 10"*, along 
with a Boussinesq simulation at 3 x 10"^, performed using 
the Kaspi et al. (2009) model. Curves give scaling predictions 
evaluated at those same values (plotted in corresponding col- 
ors). The combined scaling matches (to within a factor of 
~2) the mass-weighted mean jet speeds obtained in both the 
Boussinesq and anelastic simulations over a factor of 10^ in 
Ra*p, 10^ in Ekman number, and 10^ in Rossby number. 
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and 3 x lO""* (red). The three regimes arc clearly visible 
in the scalings and agree favorably with the simulations 
throughout the plotted range. The gradual increase in 
the slopes in Regime II near the left edge — which can be 
seen in both the scalings and the simulations, especially 
at Ekman numbers of 3 x 10~^ and 10~^ — results from 
the fact that Ra*p approaches its critical vahie (evaluated 
for the scalings using Eq. 24). Although the match is not 
perfect, the scalings not only approximately reproduce 
the asymptotic behavior suggested by Christensen but 
also exhibit the trend (suggested by Christensen's sim- 
ulations) where the transition to the asymptotic regime 
ocurs at larger Ra*p for larger E. The values of Ra*p at 
these transitions are similar in the scalings as in Chris- 
tensen's simulations to within an order of magnitude for 
most E values. Moreover, the actual values of Rossby 
number predicted by the scaling at a given Ra*p and E 
match those of Christensen's simulations to within a fac- 
tor of ~2. 

Figure 106 compares our combined scaling against sim- 
ulations performed using the Kaspi et al. (2009) model. 
Magenta squares, blue upward-pointing triangles, and 
red right-pointing triangles present the Rossby num- 
bers from anelastic simulations performed respectively at 
£■ = 1.5 X 10-4, 4 X 10-4, and 8 x lO^^ using the Jovian 
radial structure as in Kaspi et al. (2009), while the black 
squares present Rossby numbers from a Boussinesq sim- 
ulation at = 3 X IQ-^ performed with the same model. 
Scaling predictions are overplotted in the corresponding 
colors (black, magenta, blue, and red for Ekman num- 
bers of 3 X 10-^ 1.5 X 10-4, 4 X 10-4^ and 8 x 10^4^ 
respectively). In these runs, the Prandtl number equals 
10, leading to a smaller critical Raylcigh number than 
given in Eq. (24); as a result, the plotted ranges are suf- 
ficiently supercritical that the scaling predictions do not 
exhibit a gradual change in slope near the left edge of the 
plot (in contrast to the behavior in panel a). Again, the 
scalings agree well with the simulations; the discrepancy 
is everywhere less than a factor of 2. This agreement is 
encouraging considering that, together. Fig. lOo, and b 
present results spanning factors of 200,000 in Ra*p, 1000 
in Ro, and 80 in E. 

Importantly, when the modified flux Rayleigh number 
is defined in the appropriate mass- weighted way (Eq. 7), 
the same scalings explain the mass-weighted, global-mean 
Rossby numbers of both the Boussinesq and anelastic sim- 
ulations, despite the fact that the former use constant 
mean density while the latter adopt background densities 
(and thermal expansivities) that vary radially by several 
orders of magnitude from top to bottom. The Boussi- 
nesq simulations exhibit relatively little axial shear of 
the zonal wind along columns of constant .s (Christensen 
2002), while the anelastic simulations described here ex- 
hibit significant shear (Kaspi et al. 2009). The fact that 



our scaling arguments can explain both sets of runs in- 
dicates that this shear (if present) does not control the 
mass- weighted mean interior wind speed. 

Interestingly, Fig. 9 suggests that Regime I disappears 
for sufficiently small Ekman and modified flux Rayleigh 
numbers {E < 10~^ and Ra*p < 10~^ for the parame- 
ters chosen). If so, then at very small Ekman numbers, 
Regime II abuts directly against Regime III. No numer- 
ical simulations have yet been performed at these small 
parameter values, however, and how this transition oc- 
curs remains an open question. Equation (39) may need 
modification in this low-iia^ and low-E portion of the 
parameter space pending further numerical and theoret- 
ical work. 

Next, we extrapolate the scalings to Jupiter, where the 
Ekman and modified flux Rayleigh numbers are ~10-^^ 
and '^10-^^-10-^4^ respectively. Despite the uncertain- 
ties, Eq. (39) would suggest that Jupiter lies in the 
asymptotic regime. For the Jovian values of Rap and 
E, the scalings predict Rossby numbers of ^3 x 10-^- 
6 X 10-4 and the implied zonal wind speeds in the interior 
are ~0.1-1 msec-^. Although numerical simulations that 
are overforced by orders of magnitude can easily produce 
Jupiter- like jet speeds (Christensen 2001, 2002; Heimpel 
and Aurnou 2007), both Christensen's asymptotic flt and 
our mixing-length scalings therefore suggest that, given 
the weakly forced conditions of Jupiter's interior, the 
zonal-jet speeds in the interior are much slower than the 
observed cloud-level values (which reach ^150 m sec-^ for 
Jupiter and 400m sec- ^ for Neptune). Nevertheless, this 
conclusion is tentative, because our estimates of eddy vis- 
cosity are uncertain and because the asymptotic regime 
(if it exists!) is not yet obvious in numerical simulations. 
Stronger conflrmation of the existence or absence of the 
asymptotic regime in numerical models is necessary. 

7. Conclusions 

Over the past two decades, several authors have per- 
formed three-dimensional numerical simulations of low- 
viscosity convection in rapidly rotating spherical shells 
to test the hypothesis that convection in the molecular 
interiors pumps the fast zonal jet streams observed on 
Jupiter, Saturn, Uranus, and Neptune. These studies 
show that the zonal jet speeds can range over many or- 
ders of magnitude depending on the heat flux, viscosity, 
and other parameters (e.g., Christensen 2002; Kaspi et al. 
2009). Here, we presented simple theoretical arguments 
to explain the characteristic jet speeds — and their depen- 
dence on heat flux and viscosity — seen in these studies. 

Our main findings are as follows: 

• We demonstrated that the characteristic convec- 
tive velocities in the Kaspi et al. (2009) simulations 
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scale well with (agF/pCpQy^^. While such a ro- 
tating scaling has long been proposed and is con- 
sistent with rotating laboratory experiments (e.g., 
Fernando et al. 1991), this to our knowledge is the 
first demonstration of its applicability to the interi- 
ors of giant planets over a wide range of heat fluxes 
and rotation rates. 

Next, we attempted to explain the jet speeds ob- 
tained by Christensen (2002) and Kaspi et al. (2009) 
in three-dimensional numerical simulations of con- 
vection with free-slip boundaries. At weakly super- 
critical Rayleigh numbers, linear theory (e.g., Busse 
1983a; Busse and Hood 1982) and numerical simu- 
lations (e.g., Christensen 2002) show that the east- 
ward and outward convective velocity components 
are highly correlated outside the tangent cylinder. 
We showed that if the degree of correlation C be- 
tween the eastward and outward convective velocity 
components is constant (as expected at very low su- 
percriticalities) and if the jets are damped by a nu- 
merical viscosity, then the mean jet speeds should 
scale as the convective heat flux over the viscosity 
(Regime II). This scaling explains the jet speeds ob- 
tained in both Boussincsq and anclastic simulations 
at low Rayleigh numbers, where the jet speeds are 
relatively weak. 

On the other hand, at higher Rayleigh number, the 

eastward and outward convective velocity compo- 
nents become increasingly decorrelated. We showed 
that, if the fraction of convective energy release 
used to pump the jets (e) is a constant in this 
regime, and if the jets are damped by a numeri- 
cal viscosity, then the mean jet speeds should scale 
as (F jvY^"^ ■ This scaling (Regime I) explains quite 
well the mean jet speeds obtained in the simulations 
from Christensen (2002) and Kaspi et al. (2009) 
in cases when the jets dominate the total kinetic 
energy, as occurs at highly supercritical Rayleigh 
numbers. 

The relationship between the correlation coefficient 
C and the jet-pumping efficiency e shows how the 
transition between these two regimes can naturally 
occur. A nearly constant correlation coefficient im- 
plies that the jet-pumping efficiency must increase 
strongly with Ra*p. When the jet-pumping effi- 
ciency finally plateaus near its maximum value of 
1, then the correlation coefficient must begin to de- 
crease strongly with further increases in Ka*p (sec 
Eq. 28). Thus, a natural regime transition occurs 
between Regime II at low Rayleigh numbers and 
Regime I at high Rayleigh numbers. Moreover, we 
rigorously calculated C and e from the anelastic 



simulations and showed that, as expected, C de- 
creases and e increases with Rd*p. A simple, heuris- 
tic theory that transitions smoothly between con- 
stant correlation coefficient at low Ra*p and a con- 
stant jet-pumping efficiency at high Rd*p explains 
the qualitative features of the correlation coefficient 
and jet-pumping efficiency as calculated from the 
simulations. 

• Both the Boussinesq (Christensen 2002) and anelas- 
tic simulations hint at the existence of a third 
regime where, at sufficiently large heat flux and/or 
sufficiently small viscosity, the characteristic jet 
speeds become independent of the viscosity. We 
constructed a simple mixing-length scaling for the 
behavior in this regime under the assumption that 
the jet-pumping efficiency is constant and the damp- 
ing results from an eddy viscosity (rather than 
the molecular or numerical viscosity). This scal- 
ing suggests that the mean jet speed in the molecu- 
lar interior should scale with the convective heat 
flux to the 1/4 power. Given the simplicity of 
our model, this agrees favorably with the asymp- 
totic fit suggested by Christensen (2002), which 
proposes that jet speeds should scale with the con- 
vective heat flux to the 1/5 power. Both scalings 
suggest — tentatively — that the mean wind speeds 
in the molecular interior should be significantly 
weaker than the jet speeds measured at the cloud 
level. 

• A scaling that combines these three regimes can 
explain (to within a factor of ^^2) the mass-weighted 
mean jet speeds obtained in both the Boussinesq 
and anelastic simulations ranging across a factor 
of 10^ in Rd*p, 10^ in Ekman number, and 10'^ in 
Rossby number. Importantly, when the modified 
flux Rayleigh number is defined in an appropriate 
mass- weighted way, the same scalings apply to both 
Boussinesq and anelastic cases. 

It is important to emphasize that the scaling theories 
presented in Sections 4-6 apply only to the mass- weighted 

mean wind speeds. The theory docs not provide informa- 
tion on the detailed three-dimensional structure of the 
jets, such as their variation in latitude. 

Several challenges remain for the future. First, we did 
not consider the effect of the Prandtl number on the jet 
speeds, mostly because a careful characterization of the 
Ra*p- and i?-dependence of the jet speeds over a wide 
range has only been performed for a very limited number 
of Pr values (1 for Christensen (2002) and 10 for Kaspi 
et al. (2009)). However, several studies show that Pr 
does affect the jet speeds (e.g., Aubert et al. 2001; Chris- 
tensen 2002; Kaspi 2008), although one would expect this 
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dependence to disappear in the asymptotic limit. Ex- 
tending our scalings to include the Pr dependence is an 
important goal for future work. 

Second, the simulations we set out to explain — those of 
Christensen (2002) and Kaspi et al. (2009) — assume the 
convection occurs in a thick shell, for which the region 
outside the tangent cylinder dominates the total mass, 
volume, and kinetic energy. As shown by several au- 
thors, however, the dynamical mechanisms and details 
of jet pumping differ for the fluid inside and outside the 
tangent cylinder (e.g., Heimpel and Aurnou 2007). Thus, 
the scaling properties could differ for simulations in a thin 
shell (Heimpel and Aurnou 2007; Heimpel et al. 2005), for 
which the high-latitude regions inside the tangent cylin- 
der will more strongly influence mass- weighted mean flow 
properties. 

Third, we assumed that the wavemimber k — which 
represents the length scale for the variation of zonal-wind 
shear and Reynolds stress with cylindrical radius s — was 
a constant with the Rayleigh and Ekman numbers. To 
zeroth order, this seems a reasonable approximation, be- 
cause simulations suggest that the widths of the jets out- 
side the tangent cylinder are set by the adopted shell 
thickness and not by the Ekman and Rayleigh numbers. 
Nevertheless, analytic solutions in the linear limit sug- 
gest that the zonal wavelengths near the critical Rayleigh 
number should scale as E^^^ (e.g., Busse 1970; Roberts 
1968), and scaling arguments in the context of convec- 
tion with no-slip boundaries suggest that the sizes of vor- 
tex rolls scale as E^^^ (Aubert et al. 2001). Thus, the 
wavenumber representing the s- variation of the zonal-jet- 
shear and Reynolds stresses could potentially have a weak 
dependence on E. Although our assumption of constant 
k seems to work reasonably well over the simulated range, 
it would be worth relaxing this assumption in the future. 

Our scaling arguments — and the simulations they at- 
tempt to explain — assume the fluid is electrically neu- 
tral. However, the interiors of Jupiter and Saturn be- 
come metallic at pressures exceeding ~1 Mbar (Guillot 
et al. 2004), and the resulting magnetohydrodynamic ef- 
fects could influence^ th{^ jet dynamics even in the molec- 
ular region (e.g., Glatzmaier 2008; Kirk and Stevenson 
1987; Liu et al. 2008). For example, Liu et al. (2008) ar- 
gued that Ohmic dissipation would limit strong jets to re- 
gions above 96% (86%) of the radius on Jupiter (Saturn) . 
They further argued that, if the interior is isentropic and 
exhibits columnar Taylor-Proudman behavior, the winds 
would be weak not only in the electrically conducting re- 
gion but throughout the convection zone. This result is 
consistent with ours although it invokes different physics. 



In the equatorial plane, the simulations generally exhibit west- 
ward flow at the inner boundary and eastward flow at the outer 
boundary, with a relatively smooth transition in between. 



Moreover, the simulations we investigated do not include 
solar forcing or latent heating. These effects may induce 
significant internal density perturbations — and therefore 
thermal- wind shear within a few scale heights of the ob- 
servable cloud deck (e.g., Kaspi and Flierl 2007; Lian and 
Showman 2008, 2010; Schneider and Liu 2009; Wilhams 
2003). The development of convective models that in- 
clude these effects is a major goal for the future. 
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